program mkhisr
  !
  ! Listmode to Histogram Converter
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  use suty
  implicit none
  integer,parameter::MAXPATHLEN=1024
  type t_block_pair
     sequence
     integer id1
     integer id2
     integer rot
  end type t_block_pair
  type t_grouped_pair
     sequence
     integer pid
     integer gid
  end type t_grouped_pair
  integer,parameter::npid = 4
  integer,parameter::ngid = 5
  integer,parameter::md=8
  integer,parameter::mt=12
  integer,parameter::ma=12
  type(t_block_pair),allocatable::idt(:,:)
  type(t_grouped_pair),allocatable::ridt(:,:)
  integer ia,narg,len,nf,ifile
  character(MAXPATHLEN) str,fin
  character*(MAXPATHLEN),target::fout,foutp,foutd
  character*(MAXPATHLEN),pointer::s
  integer*2,target,allocatable:: hists(:,:,:,:,:,:,:)
  integer*2,pointer::histp(:,:,:,:,:,:)
  integer*2,pointer::histd(:,:,:,:,:,:)
  integer*2,pointer::hist(:,:,:,:,:,:)

  narg=command_argument_count()
  if(narg<3) then
     call usage
     stop "too few arguments"
  end if
  
  nf=0
  do ia=narg-1,1,-2
     call get_command_argument(ia,str,len)
     if(len==2) then
        select case(str(2:2))
        case("o")
           s=>fout
        case("p")
           s=>foutp
        case("d")
           s=>foutd
        case default
           call usage
           stop "*** invalid option."
        end select
        if(associated(s)) then
           call get_command_argument(ia+1,s)
           nullify(s)
        end if
        nf=ia-1
     else
        if(fout/="".or.foutp/="".or.foutd/="") then
           exit
        else
           call usage
           stop "*** please specify outfile correctly."           
        end if
     end if
  end do
 
  call setIdt
  call setRidt

  allocate(hists(mt,ma,mt,ma,npid,ngid,2))
  histp => hists(:,:,:,:,:,:,1)
  histd => hists(:,:,:,:,:,:,2)
  hists=0_2
  do ifile=1,nf
     call get_command_argument(ifile,fin,len)
     fin=trim(adjustl(fin))     
     call readList(fin,histp,histd)
  end do
  
  write(*,*) "Max Prompt: ",maxval(histp)
  write(*,*) "Max Delayed:",maxval(histd)

  if(foutp/="") call write_hist(histp,foutp)
  if(foutd/="") call write_hist(histd,foutd)

  hist=>histp
  hist=histp-histd
  where(histp<0) histp=0_2
  write(*,*) "Max True:   ",maxval(hist)

  if(fout/="") call write_hist(hist,fout)

contains
  
  subroutine write_hist(h,outfile)
    integer*2,intent(in)::h(:,:,:,:,:,:)
    character*(*),intent(in)::outfile
    integer istat
    character*128 s
    write(*,"(a,$)") "Writing "//trim(outfile)
    open(unit=1,file=outfile,access="stream",status="replace",iostat=istat)
    if(istat/=0) then
       write(*,*) "*** error during open"
    else
       write(1,iostat=istat) h
       close(1)
       if(istat/=0) then
          write(*,*) "*** error during write"
       else
          s=trim(i2a(mt))//"x"//trim(i2a(ma))
          write(*,*) "size: "//trim(s)//"x"//trim(s)//"x"//trim(i2a(npid))//"x"//trim(i2a(ngid))&
               //trim(i2a(sizeof(h(1,1,1,1,1,1))))//" = "&
               //trim(i2a(sizeof(h)))//" bytes"
       end if
    end if
  end subroutine write_hist

  subroutine readList(f,h_i,r_i)
    character*(*),intent(in)::f
    integer*2,intent(inout)::h_i(:,:,:,:,:,:)
    integer*2,intent(inout)::r_i(:,:,:,:,:,:)
    integer istat
    integer*4 dws(2)
    logical prompt
    integer t1,a1,d1
    integer t2,a2,d2
    integer nread,nacc,nvalid
    integer ip,ig

    write(*,*) "Opening file: "//trim(f)
    open(unit=1,file=f,status="old",access="stream",iostat=istat)
    stop "Failed to open the file"

    nread=0
    nvalid=0
    nacc=0
    read: do while(istat==0)
       read(1,iostat=istat) dws(1)
       if(istat/=0) exit
       nread = nread + 4

       if(btest(dws(1),31)) cycle
       if(.not.btest(dws(1),30)) cycle

       isup: do 
          read(1,iostat=istat) dws(2)
          if(istat/=0) exit
          nread = nread + 4

          if(btest(dws(2),31)) cycle read
          if(btest(dws(2),30)) then
             dws(1)=dws(2)
             cycle isup
          else 
             exit isup
          end if
          
          nvalid = nvalid + 1

          prompt=btest(dws(2),28)

          t1=ibits(dws(1),0,6) +1
          a1=ibits(dws(1),6,6) +1
          d1=ibits(dws(1),12,4)+1
          
          if(t1>mt.or.a1>ma.or.d1>md) cycle

          t2=ibits(dws(2),0,6) +1
          a2=ibits(dws(2),6,6) +1
          d2=ibits(dws(2),12,4)+1

          if(t2>mt.or.a2>ma.or.d2>md) cycle

#ifdef DEBUG
          if(green(OPT_LMF_FLIP_PMTID)) then
             if(d1/=1) d1=md-d1+2
             if(d2/=1) d2=md-d2+2
          end if
          
          if(green(OPT_LMF_SWAP_AT)) then
             call swap(a1,t1)
             call swap(a2,t2)
          end if

          if(green(OPT_LMF_FLIP_AXIAL)) then
             a1=ma-a1+1
             a2=ma-a2+1
          end if
          
          if(green(OPT_LMF_FLIP_TAXIAL)) then
             t1=mt-t1+1
             t2=mt-t2+1
          end if
#endif

          if(d1>d2) then
             call swap(a1,a2)
             call swap(t1,t2)
             call swap(d1,d2)
          end if

          ip=ridt(d1,d2)%pid
          ig=ridt(d1,d2)%gid
          if(ip==0.or.ig==0) cycle
          if(ip<0) then
             call swap(a1,a2)
             call swap(t1,t2)
             ip=-ip
          end if
          nacc = nacc +1

          if(prompt) then
             h_i(t1,a1,t2,a2,ip,ig)=h_i(t1,a1,t2,a2,ip,ig)+1_2
          else
             r_i(t1,a1,t2,a2,ip,ig)=r_i(t1,a1,t2,a2,ip,ig)+1_2           
          end if

       end do isup
    end do read
    close(1)

    write(*,*) "bytes read: "//i2a(nread)
    write(*,*) "valid/accepted = "//i2a(nvalid)//" / "//i2a(nacc)
    write(*,*) "Prompt max/sum = "//trim(i2a(maxval(h_i)))//" / "//trim(i2a(sum(h_i)))
    write(*,*) "Delayed max/sum = "//trim(i2a(maxval(r_i)))//" / "//trim(i2a(sum(r_i)))
  end subroutine readList
  
  subroutine setRidt
    integer ip,ig
    integer id1,iip
    type(t_block_pair) bp
    pointer(pbp,bp)
    if(allocated(ridt)) deallocate(ridt)
    allocate(ridt(md,md))
    ridt%pid=0
    ridt%gid=0
    do ig=1,ngid
       do ip=1,npid
          pbp=loc(idt(ip,ig))
          id1=idt(ip,ig)%id1
          iip=ip
          if(id1<0) then
             id1=-id1
             iip=-iip
          end if
          ridt(id1,bp%id2)%pid=iip
          ridt(id1,bp%id2)%gid=ig           
       end do
    end do
  end subroutine setRidt
  
  subroutine setIdt
    
    if(allocated(idt)) deallocate(idt)
    allocate(idt(npid,ngid))
    
    idt%id1=0
    idt%id2=0
    idt%rot=0
    
    ! group 1, dd=0
    idt(1,1)%id1=1; idt(2,1)%id1=2; idt(3,1)%id1=3; idt(4,1)%id1=4
    idt(1,1)%id2=5; idt(2,1)%id2=6; idt(3,1)%id2=7; idt(4,1)%id2=8
    idt(1,1)%rot=0; idt(2,1)%rot=1; idt(3,1)%rot=2; idt(4,1)%rot=3
    
    ! group 2, dd=1
    idt(1,2)%id1=1; idt(2,2)%id1=3; idt(3,2)%id1=5; idt(4,2)%id1=-2
    idt(1,2)%id2=4; idt(2,2)%id2=6; idt(3,2)%id2=8; idt(4,2)%id2=7
    idt(1,2)%rot=0; idt(2,2)%rot=2; idt(3,2)%rot=4; idt(4,2)%rot=6
    
    ! group 3 dd=-1
    idt(1,3)%id1=2; idt(2,3)%id1=4; idt(3,3)%id1=-1; idt(4,3)%id1=-3
    idt(1,3)%id2=5; idt(2,3)%id2=7; idt(3,3)%id2=6 ; idt(4,3)%id2=8  
    idt(1,3)%rot=1; idt(2,3)%rot=3; idt(3,3)%rot=5 ; idt(4,3)%rot=7  
    
    ! group 4, dd=2
    idt(1,4)%id1=1; idt(2,4)%id1=3; idt(3,4)%id1=5; idt(4,4)%id1=-1
    idt(1,4)%id2=3; idt(2,4)%id2=5; idt(3,4)%id2=7; idt(4,4)%id2=7 
    idt(1,4)%rot=0; idt(2,4)%rot=2; idt(3,4)%rot=4; idt(4,4)%rot=6 
    
    ! group 5, dd=-2
    idt(1,5)%id1=2; idt(2,5)%id1=4; idt(3,5)%id1=6; idt(4,5)%id1=-2
    idt(1,5)%id2=4; idt(2,5)%id2=6; idt(3,5)%id2=8; idt(4,5)%id2=8 
    idt(1,5)%rot=1; idt(2,5)%rot=3; idt(3,5)%rot=5; idt(4,5)%rot=7 
    
  end subroutine setIdt
  
  subroutine usage
    write(*,100)
100 format(x, "Usage: mkhis infiles",/,&
         8x,"{-o outfile_true | -p outfile_prompt | -d outfile_delayed}")
  end subroutine usage
  
end program mkhisr
